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ABSTRACT 

We derive general equations for axisymmetric Newtonian MHD and use these as the 
basis of a code for calculating equilibrium configurations of rotating magnetised neu- 
tron stars in a stationary state. We investigate the field configurations that result from 
our formalism, which include purely poloidal, purely toroidal and mixed fields. For the 
mixed-field formalism the toroidal component appears to be bounded at less than 7%. 
We calculate distortions induced both by magnetic fields and by rotation. From our 
non-linear work we are able to look at the realm of validity of perturbative work: we 
find for our results that perturbative-regime formulae for magnetic distortions agree 
to within 10% of the nonlinear results if the ellipticity is less than 0.15 or the aver- 
age field strength is less than 10 17 G. We also consider how magnetised equilibrium 
structures vary for different polytropic indices. 

Key words: stars: neutron - magnetic fields — gravitational waves 



1 INTRODUCTION 

The physics of neutron stars classes them among the most extreme objects in the known Universe: their densities, rotation 
rates and magnetic fields are all among the highest known for any astrophysical object. Typical neutron star magnetic 



fields are up to 



10 1 



G, whilst for magnetars this figure is ~ 10 ° G. Since magnetic fields induce a distortion in a 



star, a rotating magnetised neutron star could be a significant source of gravitational radiation. With the advent of second- 
generation gravitational wave detectors like Advanced LIGO, we may soon be in a position to observe neutron stars through 
their gravitational radiation signals — and hence have a new probe of the physics of these stars. 

Understanding magnetic distortions requires an understanding of the neutron star's interior field; NSs with relatively 
weak exterior fields could still have significant ellipticities if they have a much stronger field in their bulk. Here we model a 
NS as an infinitely conducting polytropic fluid and examine various kinds of magnetic field: purely poloidal, purely toroidal 
and mixed-field configurations. The numerical scheme we use is able to deal with extremely strong fields and fast rotation, so 
we are able to study the theoretical properties of very highly magnetised stars, as well as examining how well perturbative 
results hold away from the weak-field regime. 

It has long been predicted that magnetic fields will distort a fluid star ( Chandrasekhar fc Fermil 1953h . It was found that 
this distortion only becomes appreciable if the magnetic energy £ mag of the star is comparable with its gravitational energy 
W; since neutron stars have tremendous self-gravity it follows that one would only expect very strong magnetic fields to 
generate any significant distortion. This suggests that one would expect magnetars to be the most distorted NSs and hence of 
most interest to gravitational wave astronomy (with the caveat that this early work is for an incompressible fluid and so is of 
limited relevance to NSs) . For the results presented in this paper we will quantify this statement by scaling our code-generated 
results to real neutron star values. 

A number of studies of magnetically deformed stars exist. These have included work focussed on poloidal, toroidal or 
mixed fields, and boundary conditions where the fields either vanish on the surface of the star or decay at infinity. Changing 
any of these can lead to very different results, so the uncertainty we have about the geometry of NS magnetic fields translates 
into an uncertainty about how distorted neutron stars are. 

Analytic approaches have been restricted to weak fields and small deformations, as the nonlinear nature of stronger 
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magnet i c field s rapidly makes t he problem intrac t able. Early work treated deformations of incompressible fluids (see, e.g., 
Ferrarol (| 19541 ); iRobertsl (|l955l ); lOstriker fc Gunnl jl969h ). a simplifying assumption but not terribly physical for real stars. 



The first studies of compressible s tars ass umed very simp listic density distributions and magnetic fields confined within the 
star (jWoltierl Il960l ; IWentzellll96ll): later iGoossensf l|l972h treated the problem of a poloidal field matched to an external 
dipole, extending the work of Ferraro ] ( 1954 ). More re cently, work has focussed on the problem of magnetic deformations 
related specifically to neutron stars (|riaskell et~aL 20081 ). including a mixed-field case with vanishing exterior field. 

In add i tion to an alytic work , a nu mber of studies have used numerical methods to calculate magnetic distortions. 
Monaghan |l965l ) and iRoxburghl (| 19661 ) calculated field geometries and surface distortions for various polytropes, allow- 



ing for an exterior magnetic field. Their work was perturbative and so restricted to weak field s. More recently a second-order 
perturbation technique has been applied for the strong fields found in magnetars (|lokall200ll ). Other studies of highly mag- 
netised stars have solved the fully non-linear problem, to allow for more highly deformed configurations than could be accu- 
rately determined using a p erturbative approach. This was originally done for strong magnetic fields confin ed within the star 
JOstriker fc Hartwick 19681 ). by extending an earlier self-consistent field method for rapidly-rotating stars ( Qstriker fc Mark 
19681 ). For purely poloidal f ields an improve d numerical method was devised which enabled the calculation of highly distorted 



equilibrium configurations ( Miketinac 19751 ); it was found that for very strong fields the maximum density of the star could 
move away from the c entre to make the geometry of the density distribution toroidal. Solutions have also been found using a 
mixed-field formalism ( Tomimura fc Eriguchilboosl ). Finally, relativistic effe cts have been considered : fully relativistic solutions 
for purely poloidal fields (|BoTmiBt''et**al' 1 19951) and purely toroidal fields llKiuchi fc Yoshida 20081 ) and partially-relativistic 
solutions in the mixed-field case ( Kiuchi fc Kotake 20081 ; Colaiuda et al. 20081 ). In the Discussion we shall return to the role 
of boundary conditions in the mixed-field case. 

This paper is a study of the various stationary, axisymmetric equilibrium solutions for Newtonian fluid stars in perfect 
MHD. We show that the full equations of MHD reduce under these limits to two general cases: a mixed-field case ( which includes 
surely poloidal fields as a speci al case) and purely tor o idal f ields. The mixed-field formalism dates back to iGrad fc Rubin 



ljl958l ) and was recently used by Tomimura fc Eriguchil (2005) to study mixed- field stars. We are not aware of any previous 
work using the other, purely toroidal, case for Newtonian MHD. In the mixed-field case the toroidal fields vanish outside the 
star, but the poloidal fields only decay at infinity; we consider this boundary condition more realistic than the condition of 
zero exterior fields used by much of the previous work discussed above. With our formalism, we investigate the resulting field 
configurations, including the relative strengths that toroidal and poloidal fields can have, and the maximum theoretical field 
strength a fluid star can have whilst remaining in an axisymmetric stationary equilibrium state. We also look at distortions 
induced by magnetic fields, including the effect of changing the polytropic index. We examine the validity of perturbative 
results for magnetic distortions in the strong-field regime. Finally, we rescale all our code results to canonical neutron star 
values and ensure we are always comparing magnetic and rotational effects in the same physical model star. 



2 AXISYMMETRIC FORMALISM 
2.1 Governing equations 

We model a rotating magnetic neutron star by assuming that it is in a stationary state, axisymmetric with both the magnetic 
dipole axis and the spin axis aligned, and comprised of infinitely conducting material (the perfect MHD approximation) . We 
work in electromagnetic units. The equations that describe this system are the Euler equation 

_ lyp - v$ 9 + V$ r + - = 0, (1) 
p p 



together with Poisson's equation 
Ampere's law 

and the solenoidal constraint 



A$ 9 = 4ttG P , (2) 
V x B = 47rj (3) 



V • B = 0. (4) 

We close the system of equations by assuming a barotropic equation of state: 

P = P{p). (5) 

In the above equations P, p, $ 9 , <3> r , j, B, G and C are the pressure, density, gravitational potential, centrifugal potential, 
current density, magnetic field, gravitational constant and Lorentz force (£=j X B), respectively. 

Although the formalism allows for different choices of the centrifugal potential $ r and equation of state P — P(p), we 
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will work with a rigidly rotating star: 



*r=5f", (6) 

where the angular velocity flo is a constant and zu the cylindrical polar radius; and a polytropic equation of state: 

P = fcp 1+1/JV , (7) 

where k is some constant and N the polytropic index. 

The assumption of axisymmetry simplifies the equations considerably. Taking the curl of equation (JTJ leaves us with the 
requirement that 

Vx(f)=0. (8) 
Additionally, the solenoidal nature of B allows us to write it in terms of some streamfunction u, defined through the relations 

B„ = = !j£.. (9) 

zu oz zu ozu 

One may also define a solenoidal B-field by using the vector potential A, where B = V x A; we will use the (^-component A$ 
later. These two definitions are related by u — zuA^. We also define a differential operator A, by 

d 2 Id d 2 

A *= <TT--^ + #2- (10) 

ozu zu ozu oz' ! 

Using the two conditions ([5| and @, one can show that Ampere's law in axisymmetry may be rewritten as 

47tj = — X/izuBj,) x e<t, - — A»u e (11) 

zu zu 

— see section [ATI for a full derivation. 



2.2 Mixed- field formalism 



In axisymmetric perfect MHD with mixed poloidal and to roidal fields, the magnetic field and current are related through the 
Grad-Shafranov equation (see, e.g.. lGrad fc Rubin (1958) or section IA2I1 : 



iirp 



dM 



A,« + f(u 



df 



du zu 2 \ du 

where f(u) = zuB^ and M(u) is defined through VM(ii) = C/p. Combining the Grad-Shafranov equation with (|11[) yields 



(12) 



1 df dM 

J= 47^ B + ^^- (13) 
Finally we use the notation of Tomimura fc Eriguchi (2005) and Chandrasekhar fc Prendergast ( 19561 ). making the replace- 
to arrive at our final expression relating the current and field: 

j = q(u)B + zupK(u)ej, (14) 



ments a 



1 df 



and k : 



dM 
du ' 



for a mixed poloidal and toroidal field in axisymmetry. 

The two functions a(u) and k(u) govern different aspects of the magnetic field: firstly, since £ — j x B we have C — 
zupKe^, x B (from equation (|14[l ) — i.e., the Lorentz force is dependent on k, and so k governs the relative contributions of 
the magnetic and centrifugal forces to the overall distortion of the star. The role of a is less clear. From equation (|14|l we 
see that a = gives a purely toroidal current and hence poloidal field, whilst increasing a increases the size of the mixed 
toroidal-poloidal term aB (and so indirectly increases the toroidal component of the field). However, there is no limit in which 
the field is purely toroidal in this formalism. We can thus only expect a to have some indirect connection with the relative 
strengths of t he poloidal and toroida l comp onents of the magnetic field. 

Following Tomimura fc Eriguc hi (2005), we choose the functional forms of a(u) and k(u) as: 



k(u) 



ko 



const. 



a(u) 



Oj{u U m ax^ if U ^> Umax 







if 11 ^ Umax 



(15) 



(16) 



where a is chosen to ensure there is no exterior current, £ is some constant and u max is the maximum surface value attained 
by the streamfunction u. Next we combine the definitions a = ^f- and f(u) = zuB^ to see that 

/u 
q(m') du = zuB^ (17) 

— i.e., we must enforce the continuity of f a(u) du to ensure the continuity of B^. We therefore choose the lower limit of the 
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integral of a so that 

/u 
q(m') du 



(^ ^mai) if If ^> Umax , (18) 

if M < M maI . 



For our chosen functional forms of q(u) and k(u) we see that for a specific solution we need to choose three constants: 
a and kq. We will later drop the zer o subscript, with the understanding that k always refers to a constant unless otherwise 



stated. iTomimura fc Eriguchil (|2005l ) set £ = 1, but we have found that a smaller value of £ allows for a slightly stronger 
toroidal-field component; accordingly, we set £ = 0.1 throughout this paper, except in comparing our results with previous 
work (subsection 13. 4|l . 

For the purposes of numerics we seek integral equations; the integral form of JT} is 

H = C -<S> g + $,.+ f \(u')du', (19) 
Jo 

where C is an integration constant and 



(20) 



(1 



is the enthalpy. 

The integral form of Poisson's equation ((2)1 is: 



<M r ) = ~ G I lT^TT (21) 



r — r' 



Finally, using the current relation (|14|l an integral equation for the magnetic field may be found (see lTomimura fc Eriguch 

( 20051 ) for details): 



. , . f 4ty r- 2 - In a(u) au + Kpzu , . , 

A Jr) sin cj>= ^ ^- — sin ^' dr. (22) 

J l r - r 

With the three equations l)19p. (|21|l and (| 22f) it is possible to calculate stationary configurations of magnetised rotating 
stars (together with the specified constants a and At). 

2.3 Toroidal-field formalism 

For a purely toroidal field we have B = B^e^. In this case, the A*« term disappears from equation (|ll|l and the Lorentz force 
reduces to the form C = B^j po i x e^. Comparing these two expressions, one can show that B^ is related to 7 = pvj 1 through 
some arbitrary function h (which must vanish outside the star): 

B^ = ifc( 7 ) (23) 

VJ 

— see section E4l for details. The magnetic potential for this case is: 

1 f™ 2 hdh , 

M = -— / -— dry, 24 
4tt J 7 d7 

where VAf = C/p as before. 

For simplicity we choose h{pzu' 2 ) — Xpvj 2 where A is a constant we specify for each code run. With this choice of h we 
then have B^, = \pzu and M = — X 2 pzu 2 , so that the first integral of the Euler equation becomes 

H = C ~$+\Q. 2 vj 2 - ±\ 2 pzu 2 . (25) 

This equation, together with Poisson's equation, is sufficient to find numerical solutions; the toroidal-field case is thus simpler 
than the mixed-field formalism, which also had an extra equation for the magnetic field. 

2.4 Restrictions on the magnetic functions 

In the appendix (outlined in the above sections) we show that for axisymmetric perfect MHD in a fluid, the equations reduce 
to a mixed-field case (with two magnetic functions a(u) and k(u)) and a purely toroidal case (with a magnetic function h(-y)). 
Although the magnetic functions appear to be arbitrary, there are a number of restrictions on their functional forms, on either 
physical grounds or because they result in trivial solutions. 

The functions a(u) and ^1(7) (where u = vaA,/, and 7 = pvj 2 as before) govern the toroidal fields in the two cases, 
and so both must necessarily vanish o utside the star to avoid havin g exterior currents. Since the streamfunction u does not 



vanish at the star's surface, we follow ITomimura fc Eriguchil ((2005) in defining u m ax to be the maximum surface value of u 



and then choose a(u) to be a power of (u — u max ), which does vanish outside the star. There does not appear to be any 
other functional form for a which vanishes outside the star and is dependent only on it, so we conclude that (|16p is the only 
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acceptable choice for a(u). The functional form of h, similarly, appears restricted. To vanish outside the star /i(7) cannot 
contain a constant piece, so let us consider a functional form of /i(7) = \*y x where A and \ are constants. However, if \ < \ 
then = A7 x m _1 = \p x zu 2x ~ 1 will diverge at the origin, so we discard these choices. Additionally, we find that if x > 1 is 
chosen, then the field iterates to zero in our numerical scheme, leading us to choose h(-y) = A7. 

Finally, the function k(u) is theoretically allowed to depend on the streamfunction it, but if it is chosen as anything other 
than a constant then, as for /i(7), we find that the configuration iterates to a zero-field solution. This may be a limitation of 
our numerical scheme rather than a physical restriction, but in either case our HSCF-scheme solutions are limited to those 
with k being equal to some constant. 

We conclude from this that, in fact, the choices made for our functional forms are not specialised ones and (at least within 
our scheme) do not result in the exclusion of physically valid solutions. Rather, we believe that our results are quite generic 
to perfectly conducting polytropes in axisymmetry. 



3 NUMERICS AND CALCULATING VARIOUS QUANTITIES 
3.1 Numerical scheme 

Ou r code uses the Hachi s u self -consistent field (HSCF) method (jHachisu (1986). extended to magnetised configurations 
by iTomimura fc Eriguchil (|2005l )) to iteratively find a stationary solution to the hydromagnetic equilibrium equation ([T]). 
Specifically, one specifies a polytropic index N, magnetic functions a(u) and re(u) and the ratio of polar to equatorial radii 
r p /r eq , and the code determines the angular velocity, density distribution and other quantities consistent with the user's input 
parameters. 

The iterative steps for our extended (mixed-field) HSCF scheme are: 



1. Make an initial guess of p=const.; 

2. Find <£> s from Poisson's equation (|21[) ; 

3. Guess ^4^=const.; 

4. Find an improved form of A$ from equation (|22|l and the earlier guesses for p and A$ (this is the iterative step for A$)\ 

5. Find fig and C from the boundary condition that the enthalpy must vanish at the surface of the star; this requires the 
potentials $ 9 and A^ found earlier and a user-specified axis ratio r p /r eq ; 

6. We now know all right-hand side terms in (|19[l ; use the equation to determine the enthalpy at all points in the star; 

7. Find the new (improved) estimate for the density distribution using p new (r) = ^ ) where N is the polytropic index 
and H m ax the maximum value of enthalpy attained in the star; 

8. As the iterative step, return to step 1 but US6 O — Onew instead of the earlier density distribution (p=const. for the first 
cycle). At step 3 in the new cycle, use the 'new' form of A^ calculated in step 4 of the previous cycle. 

This sequence of steps is repeated until the code has achieved satisfactory convergence in H max , Q$ and C. The toroidal-field 
scheme is similar to the one above except that the magnetic field is directly related to the density by = Xpzu for pure 
toroidal fields. For this reason there is no separate iteration in the magnetic field and steps 3 and 4 are no longer needed. 
The magnetic field only enters in step 6, where the enthalpy H is found from the pure-toroidal equation (|25|l instead of the 
mixed- field version (1191 1, 



3.2 Magnetic energy and field strength 

We will wish to calculate the magnetic energy £ mag of the star in the code, to compare different configurations and also to 
calculate a virial test (see section 13.31 and figure [1} . The familiar definition of £ m ag is 

f B 2 

£mag = / — dr, (26) 
all space 

but this is not suited to numerical evaluation, since the integrand only decays at infinite distance; our numerical integration 
is over a finite radius and so this definition would introduce truncation error. However a physically equivalent definition for 
£mag, more useful here, is 

£ mag = J r ■ C dr (27) 

all space 

— since C has compact support (through its p dependence) the above integrand will vanish outside the star. 
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(a) 



vc 



(b) 



VC 



MP 



10000 100000 

MP 



(C) 



VC 



MP 

Figure 1. Convergence tests for: (a) a purely poloidal field, no rotation, axis ratio of 0.2; (b) mixed field with 2% toroidal field, no 
rotation, axis ratio of 0.2; (c) purely toroidal field, no rotation, axis ratio of 1.05. Here VC is the virial test result and MP the number 
of mesh points. Since VC decreases as MP increases, we see that the code is convergent. 



For a measure of the magnetic field strength of the star, we define a volume-averaged magnetic field B through 

8ty£ ma a 



B 2 



hi' 

all space 



V 



(28) 



3.3 Virial test 

We may use the scalar virial theorem (see, e.g., Shapiro fc Teukolskv ( 198St )) as a test of convergence for the code. For a 



rotating magnetised self-gravitating fluid the virial theorem states that 



IH-2T + S 
2dt* ~ ZI +t " 



3U „ r 



(29) 



where / is the moment of inertia about the rotation axis and T,£ mag ,U and W are the kinetic, magnetic, internal and 
gravitational energies, respectively. For our stationary star I has no time variation so the first term is zero. Given this, we 
expect the various energies for our star to satisfy 



VI 

2T + £ mag + — + W = 0. 
N 



(30) 



Calculating the quantity on the left-hand side of the above equation tells us the absolute deviation from zero, but we need 
to know the relative error. A value of 2T + £ mag + 3U/N + W = 10~ 5 would appear to indicate acceptable accuracy, but if 
the individual energies are of order 10 -4 then the relative error is unacceptable: around 10%. For this reason we normalise by 
dividing through by W and define our virial test result VC as 

\2T + £ mag + 3U/N + W\ 



VC: 



\w\ 



(31) 



— the smaller the value of VC, the greater the code's accuracy. We use VC in our convergence testing, figure [T] In the figure 
we see that as grid resolution is increased the virial test result decreases; in particular, since the gradient of each plot is 
approximately —1 we conclude that the code is first-order convergent. 
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Table 1. Dimcnsionlcss quantities for a sequence of stars with TV = 1.5, k = 0.4, a = 200 and f = 1. 





£ mag /\W\ 


u/\w\ 


T/\W\ 


\W\ 


(r- 


M 


vc 


0.588 


0.144 


0.284 


1.21e-03 


4.81e-02 


5.14e-04 


0.831 


2.97e-05 


0.55 


0.151 


0.276 


l.lle-02 


4.59e-02 


4.53e-03 


0.811 


3.10e-05 


0.50 


0.165 


0.264 


2.11e-02 


4.32e-02 


8.01e-03 


0.787 


3.33c-05 


0.45 


0.189 


0.255 


2.27e-02 


4.01c-02 


7.72e-03 


0.763 


3.63e-05 


0.40 


0.222 


0.252 


1.19e-02 


3.58e-02 


3.45e-03 


0.729 


4.02e-05 


0.371 


0.242 


0.252 


1.10e-03 


3.31e-02 


2.89e-04 


0.705 


4.32e-05 



3.4 Comparison with previous work 



Tomimura fc Eriguchi Their results are nondimen- 



As a confirmation of our results, we compare with table 4 from 
sionalised by dividing by appropriate powers of p m ax, r eq and and these dimensionless quantities are denoted by a hat; 
for example 

& = A " 2 ■ (32) 

For comparison with their results we must also use ( = 1 instead of £ = 0.1 as the exponent in the functional form of a from 
(|16[) . Taking this into account we find that for a N = 1.5 polytrope, with k — 0.4 and a = 200, we have the sequence of 
configurations given in table [T] 

Our highest and lowest axis ratios (0.588 and 0.371) differ slightly from those of 
0.589 and 0.372), so we cannot make a direct comparison for these values. However for the other four axis ratios our values 
agree to within ~ 8% for Cl 2 and T/|W|, and to within ~ 1% in all other quantities. We also show our results for the virial 
test, showing that all our results have relative errors of ~ 10 . 



Tomimura &i Eriguchil (|2005l ) (who have 



3.5 Toroidal and poloidal energies for the mixed case 

The code variables k and a are related to the ratio of toroidal to poloidal field strength, but in a very nontrivial manner. To 
get a more intuitive, physical, measure of their respective strengths we would like to know the part of the magnetic energy 
contained in the poloidal and toroidal fields, £ po i and £tor, respectively. 
Since the total magnetic energy is given by 

£ mag = ^~ J B ■ B dr = ^- J (Bl + Bl + B 2 ) dr, (33) 

we define the poloidal energy by 

£ poi = f {Bl + B 2 Z ) dv = i J {Bl + B 2 Z ) drdp (34) 

and the toroidal energy by 

£ tor = ^ J Bldr=lJ^ J Bl drdp (35) 

where the integration here is over spherical polars r and /x = cos 9. Note that since B$ — outside the star, the toroidal-energy 
integral only needs to be evaluated over the stellar interior. For our mixed-field configurations we use Stor/Smag as a measure 
of the proportion of toroidal field. 



3.6 Ellipticity 

For the code we specify the axis ratio r p /r eq , which measures the distortion of the star's surface. For a measure of the 
distortion of the whole mass distribution of the star we define an ellipticity e through the (unreduced) quadrupole moments 
at the equator I eq and the poles I p : 

e= Ie \~ Ip . (36) 



3.7 Constructing physical sequences of stars 

To make a meaningful study of a group of different equilibrium configurations requires ensuring that we are always comparing 
the effects of magnetic fields and rotation in the same physical star: we do this by ensuring that we work with sequences of 
constant (physical) mass and the same equation of state — i.e. both the same polytropic index TV and polytropic constant k. 
We note that other intuitively sensible choices, for example just fixing the equatorial radius or central density, would mean 
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Figure 2. Contours of magnetic field strength. Plots (a), (b), (c) are (respectively) the poloidal, toroidal and total magnetic field strength 
for our canonical NS with a mixed-field configuration consisting of 3.0% toroidal field. Plot (d) is for magnetic field strength in a purely 
toroidal-field star. Note how the toroidal field in this case is much more extensive than in the mixed- field plot (b). In plots (a) and (c) the 
maximum field strength is 5.5 X 10 17 G (at the origin) and the contour separation is 5.5 X 10 16 G. For plots (b) and (d) the maximum 
field occurs in the centre of the torus bounding the toroidal field; the maximum values are 2.6 X 10 17 G and 2.8 X 10 17 G, with contour 
separations of 2.9 X 10 16 G for both plots. The bold red line in each plot represents the star's surface; the values on the axes show the 
nondimensional radius r/r eq (where r is the physical radius). 



comparing stars of either different mass or different equation of state. This would make quantifying the effects of magnetic 
fields and rotation more difficult. 

We fix our neutron star mass to the generic value of M — 1.4M@ = 2.8 X 10 33 g. For the equation of state we work with 
N — 1 polytropes and fix k by requiring that the radius R of the equivalent unmagnetised nonrotating (and hence spherical) 
star is 10 km. We will term this star the 'background' star, with the understanding that this refers to a configuration without 
magnetic fields or rotation, rather than having any perturbation theory connotations. Using the (N = 1) polytropic relation 
R = sJnk/2G we see that this gives a polytropic constant of 4.25 X 10 g cm°s . For the rest of this paper, when we quote 
physical parameters they will be for our 'canonical neutron star' with M = IAMq and k = 4.25 x 10 4 g _1 cm s s -2 . 

For the plots where we have used different polytropic indices the redimensionalising is less strai ghtforward, as the rela tion 
between spherical radius R and k also includes powers of the 'background' central density p c (see IChandrasekhar (|l939h for 
the required polytropic relations). For these cases we again fix the mass at IAMq and fix the background central density at 
p c = 2.19 x 10 15 g cm -3 — the same value as for the background N = 1 star discussed above. This then fixes R and k. 



4 MAGNETIC FIELD CONFIGURATIONS 

With the formalism described above, we are able to examine the field configurations generated in axisymmetric perfectly 
conducting polytropes. Since neutron star matter is thought to have high conductivity and be roughly approximated by an 
N = 1 polytrope, the field structures shown here should have some similarity to those in real neutron stars — although the 
field strengths here are considerably higher than those that have been observed so far. The plots in this section show contours 
of the magnetic field strength given by |B| = a/B ■ B, and of the poloidal and toroidal components, |B po ;| = \J B'^ + B'i 
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2.5 




Figure 3. Contours of the streamf unction u for a purely poloidal-ficld star; these contours are parallel to magnetic field lines and so 
represent the direction of the field. The surface of the star is the bold red line. Field lines for the toroidal component of a mixed-field 
star, or for purely toroidal stars, would go into the page and hence are not plotted here. 
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Figure 4. Left: the ratio of polar field to volume-averaged field, B v /B, as a function of the polytropic index N. The plot is for purely 
poloidal fields in non-rotating stars, all with an axis ratio of 0.996. Note that if the field was purely toroidal then this ratio would be 
zero, regardless of N . 



and |Bt or | = \B$\. All of the magnetic-field results presented here (and discussed in this section) are for non-rotating N = 1 
polytropes, unless otherwise st ated. In addition, we have concentrated on mixed-fie ld configurations h ere, since there are strong 
indic ations from both theory (jlVIarkev fc Tavlerlll97ol ; lTavlerlll97ol ; IWright|[l97il ) and simulations (jBraithwaite fc Nordlundl 
200rj ) that both purely poloidal and purely toroidal fields are generically unstable. 



In figure [2] we plot the poloidal (plot (a)) and toroidal (plot (b)) components of a mixed-field star and the total field of 
the configuration (plot (c)). We also plot the field structure of a star generated from our purely toroidal-field formalism for 
comparison (plot (d)). We see that the poloidal field pervades most of the interior of the mixed-field star, as well as extending 
outside it. This component of the field is highest in the centre and onl y goes to zero in a, small region at the edge of the star 
(seen as the pair of semicircular contours on the equator at x ~ 0.8); Markev fc Tavler ( 197ot ) call this zero-field point the 
'magnetic axis'. By contrast the toroidal field is wholly contained within this small region where the poloidal field vanishes; 
this region is dictated by the functional form of a (it) that we use. All configurations shown in this section are non-rotating, 
but rotation does not greatly affect the nature of the magnetic field. 

Comparing plots (b) and (d) in figure [2] we see that, although the maximum field strengths and contours are of similar 
magnitude in the two cases, the field in the pure-toroidal case extends over a far larger region of the star than the toroidal 
part of the mixed-field configuration. 

Figure [2] shows the magnitude of the magnetic field at a particular point; in figure [3] we show the direction of a typical 
poloidal field by plotting contours of the streamfunction it. These contours are parallel to magnetic field lines, from section 
of the appendix. Since a purely toroidal field has direction vector e^, the field lines would go into the page in the x — z 
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Table 2. Comparing parameters related to the influence of the toroidal component in a mixed-field star with axis ratio 0.9. 



Qi &tor j &i 



mag 



£mag/\W\ 



e 



B p /B 




10 
20 
30 
40 
50 



0.00 
9.87e-03 
3.02e-02 
3.96e-02 
4.05e-02 
3.86e-02 



2.43e-02 
2.55e-02 
2.82e-02 
2.93e-02 
2.92e-02 
2.88e-02 



0.216 
0.216 
0.213 
0.204 
0.196 
0.191 



0.580 
0.554 
0.504 
0.484 
0.488 
0.495 



plane we employ here (they would form concentric circles in the x — y plane) . Mixed-field lines lie in neither plane so we have 
not shown them here. 

Lastly in this section, figure [4] shows the dependence of the ratio B p /B on the polytropic index N; we find that there is 
an approximately linear relationship between the two, and for all polytropic indices B p /B is of the same order of magnitude. 
For N = 1, Bp/B « 0.5, suggesting that neutron stars (approximated as iV = 1 polytropes) with purely poloidal fields are 
likely to have a B around double the polar field B p . 

4.1 The relationship between a and £tor/£mag 

As mentioned earlier, we can increase the proportion of toroidal field in the mixed-field configurations only indirectly, by 
varying the code parameter a from equation (|16|l . In table [2] we show the effect of changing this parameter, for a non-rotating 
star with axis ratio r p /r eq = 0.9. One would expect that increasing a would increase the toroidal portion of the field, which in 
turn would lead to a decrease in oblateness (since toroidal fields induce prolate distortions) ; one would also expect a reduction 
in the ratio B p /B (since more of the field is toroidal and hence does not extend outside the star). Looking at the table, we see 
all of these effects do occur as the value of a is increased, up until the a — 40 configuration. At this point the larger value of a 
is no longer reflected in stronger toroidal-field effects. In all cases changing a does not strongly affect the value of £ m ag/\W\, 
confirming our expectation that it is the variation in the toroidal component which affects ellipticity and B p /B, rather than 
simply a reduction in £ m ag/\ W\. Finally, we note that even for the highest values of a, the relative contribution of the toroidal 
portion of the field is very small — only 4% of the total here. We shall see later that this is a generic feature of our formalism 
together with our boundary condition, where poloidal fields extend outside the star but toroidal ones vanish at the surface. 



5 MAGNETIC AND ROTATIONAL DISTORTIONS 

Having looked at field configurations, we now turn to the distortions these fields produce in the star's mass distribution. For 
purely poloidal fields we confirm previous work that these fields induce an oblate distortion; the surface shapes of such stars are 
thus similar to those of rotationally distorted stars. However, the interior density distributions are very different: centrifugal 
forces tend to leave a smaller high-density central region, whilst the Lorentz force acts to pull the point of maximum density 
away from the centre into a maximum-density ring. In the extreme limit where the ratio r p /r eq — > 0, the star actually becomes 
a torus (figure [5|. For mixed fields, the effect of increasing the toroidal component is similar to the effect of adding rotation: 
it tends to push the maximum density region back to the centre — see figure [6] Note that both the mixed-field stars shown 
are oblate though, due to the dominance of the poloidal component; stronger toroidal fields tend to make stars prolate, but 
our formalism and boundary condition seem to generate mixed fields with weak toroidal components only (the 5.5%-toroidal 
field of figure [6] plot (c) is relatively strongly toroidal, within this context). 

For weak fields and small distortions, perturbation theory results suggest that the ellipticity of a star should depend 
linearly on B . With our non-linear code we are able to check this, and see how well the perturbative result holds as field 
strengths are increased; this is plotted for both poloidal and toroidal fields in figure [7] The results depart slowly from the 
linear regime to begin with, but in the poloidal-field case the field strength required reaches a peak and then decreases again, 
for increased ellipticity. This peak seems to correspond to roughly the point at which the maximum density is pulled out into 
a ring, making the star's density distribution toroidally-shaped. We speculate that for very low axis ratios (i.e. very strong 
fields), this toroidally-shaped density is a more stable, lower-energy state than one where the maximum density remains at 
the centre. 

Purely toroidal fields give prolate density distributions, although we find that the surface shape remains virtually spherical 
even for large ellipticities (i.e. strong fields). Because rotation gives rise to oblateness in stars, it opposes the effect of a toroidal 
field in a star, and the two effects can balance to give a rotating magnetised star with zero overall ellipticity. Note that in this 
case the stars will have oblate surface shapes but a spherical density distribution — see figure [8] 

Next we look at the effect of magnetic fields on the Keplerian velocity Qk — see figure [9] We find that whilst increasing 
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Figure 5. Nonrotating N = 1 polytropcs distorted by the effect of a poloidal magnetic field, with axis ratios of 0.7, 0.5, 0.2, 0.0 (from 
(a) to (d)). For increasing distortion the maximum density moves away from the centre and the density distribution becomes toroidal 
(in the sense that the maximum density moves away from the centre of the star). As before, the numbers on the axes are dimensionless, 
but for our canonical NS r eq = 11.4, 12.9, 16.2, 17.0 km for plots (a)-(d). 

the field strength causes a slight decrease in the velocity needed to cause mass shedding, this effect only becomes noticeable 
for very strong fields. It seems, therefore, that magnetic fields are unlikely to affect the stability of a star in this manner. 

We have generally presented results for an N — 1 polytrope, as this is regarded as a reasonable approximation to a 
neutron star. For our final two figures, however, we briefly investigate the effect of varying the polytropic index TV, whilst 
maintaining a mass of 1.4Mq and central density of 2.19 x 10 15 g cm" 3 in the corresponding unmagnetised 'background' 
polytropic star. In figure 1101 we plot four stars with the same surface distortion r p /r eq = 0.5 but different N. We see that 
when N is low the density contours are all close to the edge of the star, with a large (slightly off-centre) high-density region; 
in the limiting case N = the star is an incompressible, uniform density configuration, so all contour lines coincide with the 
star's surface. For higher values of N the high-density region becomes smaller and the low-density outer region becomes larger. 
We note that the N = 2 polytrope shown cannot be a neutron star model, however, as its maximum density of 1.79 x 10 14 g 
cm" 3 is lower than the density of heavy nuclei, po = 2.4 x 10 14 g cm -3 . 

Finally, in figure [Til we look at non-rotating stars magnetised by a purely poloidal field, with an axis ratio of 0.95. We 
plot the dependence of the field strength on polytropic index TV, finding that as N is increased a weaker field is required to 
support the same surface distortion. 

6 DISCUSSION 

To understand how strong magnetic distortions may be in highly magnetised objects like magnetars, realistic models are 
needed to study the field structure of these stars. The formalism we use in this work comes directly from the assumptions of 
axisymmetry and perfect conductivity, together with a boundary condition that the poloidal part of the field should decay 
at infinity rather than vanishing at the star's surface; we anticipate that these conditions provide a reasonable model of a 
neutron star's field. 

The general formalism of axisymmetric MHD reduces to a mixed-field case and a purely toroidal-field case, with two 
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(a) (b) 




Figure 6. Density contours in an N = 1 polytropic star with axis ratio of 0.6, with different sources of distortion. Plots (a), (b) and 
(c) are nonrotating configurations with, respectively: purely poloidal field, mixed-field with 3.4% toroidal field, mixed-field with 5.5% 
toroidal field. Plot (d) is for a purely rotationally-distorted star with no magnetic field. All stars have the canonical mass of 1.4Mq, with 
equatorial radii of 12.1, 12.5, 13.2, 14.4 km for stars (a)-(d), respectively. We note that whilst a purely poloidal field tends to push the 
maximum density away from the centre, both toroidal field components and rotation have the effect of increasing the equatorial radius 
and making the star more diffuse. 




distortion 



Figure 7. Left: a graph showing how (poloidal) magnetic distortions vary with the field strength. 1 — r p /r eq is the surface distortion, 
whilst e represents the distortion of the density distribution, as defined in equation H36I I. Note that the required field strength peaks for 
1 — r p /r eq ~ 0.6 or e ~ 0.8 and then drops slightly for more extreme distortions. For small distortions we see that there is a roughly 
quadratic dependence on the field strength. Right: toroidal-field distortions versus B 2 . In this case we only use t to gauge the level of 
distortion, as the surface shapes remain nearly spherical. 
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Figure 8. Two stars with toroidal magnetic fields. The left-hand configuration is a non-rotating star (and hence has a prolate density 
distribution), whilst the right-hand one is the same physical star but with rotation added, with an oblate surface shape but an overall 
cllipticity of zero. The average field strength in both cases is B = 2.4 X 10 17 G. 




B 



10 17 G 

Figure 9. The dependence of Keplerian velocity Qk on magnetic field strength B, for stars with purely poloidal fields. Note that an 
appreciable decrease in Qk only occurs for very strong fields. 



(mathematically) arbitrary functions in the former case and a(u)) and one in the latter (/i(7)). Despite the apparent 

freedom in choosing these functions, we found that on physical grounds only one functional form was satisfactory for each 
one; see section 12.41 We conclude that the equations we have numerically solved in this work are in fact quite general and 
that we have not excluded physically valid branches of solutions with our choices. 

Perturbative calculations in the weak- field regime have found that e depends linearly on B 2 . With the use of our nonlinear 
code we are able to investigate how well this approximation holds for larger fields and ellipticities. We can see graphically 
that the first few points from both plots in figure [7] lie in fairly straight lines and hence we deduce the relations 

^^X10-( T ^)^2X10-( T ^) 2 (37) 

for the purely poloidal case (the above relation also uses B p /B ~ 0.5 from figure |4}, and 

e tor ~ -3 x 10- 4 ( -JL-) (38) 



10 16 G, 

for the purely toroidal case; where in both cases we have used a star of mass 1.4Mq whose radius would be 10 km if 
unmagnetised. By comparing these extrapolated linear-regime formulae with our non-linear code results, we can explore how 
well perturbative results are likely to hold in a strong-field regime. We find that the linear regime given by (|37p and (|38[l 
differs by less than 10% of the actual non-linear code result (shown in figure^ provided that B < 1.5 x 10 17 G, or equivalently 
e < 0.15. Alternatively, if we allow the linear relation to differ by up to 30% from the nonlinear result, we may use the linear 
relation as an 'acceptable' approximation for B < 3 x 10 17 G or e < 0.35 (i.e. it holds for the entire range of ellipticities we 
can plot in the toroidal-field case). 

This suggests that for all known neutron star field strengths, e is likely to be linearly dependent on B 2 , to a good 
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Figure 10. Non-rotating configurations, all with a purely poloidal field and an axis ratio of r p /r eq = 0.5. Plots (a) to (d) are for 
N = 0.5,1,1.5,2 polytropes, respectively; the corresponding field strengths arc B = 7.62,4.31,2.98,1.13 x 10 17 G, the maximum 
densities are 1.67, 1.14, 0.623, 0.179 X 10 15 g cm -3 and the equatorial radii are r eq = 10.2, 12.9, 17.6, 29.6 km, respectively. 
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Figure 11. Right: the poloidal field strength B required to induce a surface distortion of r p /r eq = 0.95, plotted for various polytropic 
indices. We see that the required field is weaker for higher- TV polytropes. 



approximation. Hence perturbation theory could provide accurate predictions of NS distortions, provided the neutron star 

model used is also a close approximation to real NS physics. 

We are also able to compare our linear-regime formulae with the analytic work of lHaskell et al. (2008), who also treated 
pure poloidal fields extending outside the star and pure toroidal fields vanishing at the stellar surface (as for our work). For 
the same mass, radius and polytropic index their formulae give: 
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where B B is the surface magnetic field strength, which was assumed constant for the calculation of Irlaskell et al, (2008): we 
do not have a constant surface field so have compared with their work using the value of |B| at the pole instead. Since their 
field geometries are clearly not identical to ours, and since we had to extrapolate to obtain our formulae, we would not expect 
precise agreement. Nonetheless, we feel that the similarities show that our work makes sensible contact with perturbative 
calculations. 

From figure [7j beginning at an unmagnetised spherical star, we find that in both the poloidal and toroidal cases the 
magnetic field strength required increases for larger distortions, initially; as would be expected from perturbative work. 
However, in the purely poloidal case the field strength then peaks at t ~ 0.8, dropping slightly as e is increased further. 
Around the same point the density distribution becomes toroidal in nature — that is, the point of maximum density moves 
away from the centre and a high-density torus forms; this leads us to speculate that at e ~ 0.8 it becomes energetically 
favourable for the density to change from a spheroidal profile (as seen in the weaker-field stars, e.g. figure [5j plot (a)) to 
a toroidal one (e.g. figure [SJ plot (d)). It is clear that if the magnetic field in a star is increased beyond the peak value of 
~ 5 x fO 17 G shown in the left-hand plot of figure [7] then one of our initial assumptions must be violated. Since we cannot 
investigate the possibilities with our current code, we conclude that a hypothetical star with a field of B > 6 x fO 17 may 
either have no equilibrium solution (in which case it may lose magnetic energy until it is in equilibrium), or that there may 
be a new triaxial branch of super-magnetised solutions bifurcating from the biaxial curve at e ~ 0.8. 

We do not find a similar peaking of the field strength in the purely toroidal case, however. In this case the largest 
ellipticities we are able to calculate are around e ~ 0.35. Whilst this particular value may represent a limitation of our 
numerical scheme, we suggest that a limited range of ellipticities is a consequence of the formalism for toroidal fields in 
axisymmetry, where B is directly linked to the density p; in the mixed-field case we have a separate equation to iteratively 
solve for the magnetic field. Thus restrictions on the field geometry may restrict the size of permissible ellipticities. 

Of course, whilst the 'peak field strength' we discuss here is a theoretical upper bound on NS fields, there are probably 
other physical effects that place a lower bound than ~ 5 x 10 17 G on the maximum field. Certainly, if magnetar surface fields 
are ~ 10 G one would not expect their volume-averag ed fields to exceed ~ 10 16 G. 

We have argued that the equations we solve in this paper lead to quite general solutions for axisymmetric stars. However, 
we find that although it is possible to find solutions with purely poloidal or purely toroidal fields, the range of mixed-field 
solutions is very limited. Using £tor/£mag as a gauge of the strength of the toroidal component in a mixed-field star, we find 
that for all our stars ^ £tor/£mag < 0.07. The other extreme is of course £t or /£mag = 1 for purely toroidal fields. This 
means that although the toroidal component does have some influence in a mixed-field star (see table at the end of section 
2J, it is dominated by the effect of the poloidal field. In particular all our mixed-field stars have oblate density distributions. 

Our mixed-field stars have the boun dary condition that t he toroidal component vanishes at the surface, whilst the poloidal 
piece only decays at infinity. By contrast, Haskell et al. ( 20081 ) considered the problem of mixed-field stars where the total field 
vanished at the surface. This results in an eigenvalue problem, with all (discrete) solutions having prolate density distributions. 
Since the chief difference between our work appears to be the choi c e of b oundary condition, we speculate that our boundary 
condition favours poloidal distort ions, whilst that of Haskell et al. I 2008h favours the toroidal component. 



The numerical simulations of Braithwaitei 1 20081 ) suggest that a stable magnetic field will have 0.20 < £tor/£mag 0.95. If 
this result is directly applicable to our work then it would imply that none of the solutions that exist within our axisymmetric 
formalism are stable. However, for numerical reasons these simulations use a magnetic diffusivity t erm which is zero within 
the st ar and increases through a transition region to a high, constant value in the exterior (see Braithwaite fc Nordlundl 
( 20061 ) for details). We suggest that this transition region may favour the toroidal component of a mixed-field star; it would 
be interesting to see if a similar stability result emerges from simluations using a boundary condition more similar to ours. 

Although we regard our boundary condition as the most natural for a mixed-field fluid with infinite conductivity, neutron 
stars are not perfect conductors. In moving from the superfluid interior to the crust and magnetosph ere, it is clear that the 
resistivity of the medium increases and hence the boundary condition should be adapted to reflect this. Colaiuda et al. (2008) 
noted this and attempted to mimic more 'natural' boundary behaviour by allowing the poloidal part of the field to extend 
outside the star (as for our field), but matching the toroidal part to a surface current rather than forcing it to vanish at the 
surface. 

With no clear idea about the nature of currents on the surface of neutron stars, we suggest that it may be easier to 
neglect their effects, so that the toroidal-field component vanishes at the surface. Incorporating the effects of resistivity in 
the outer regions of the neutron star would then involve adapting the boundary condition for the poloidal component; this 
could resemble a su rface treatment som ewhere between ours (where the poloidal field is unaffected by passing through the 
surface) and that of lHaskell et al.l (|200S ) (where the poloid al field decays at the surface). Since our boundary condition gives 
a poloidal-dominated field and that of Haskell et al. ( 2008h gives a toroidal-dominated field, we suggest that the inclusion of 



resistivity would result in configurations where neither component is universally dominant. In particular, we would not expect 
magnetic distortions in real, mixed-field, neutron stars to be universally oblate or prolate. We conclude that future, more 
realistic, models of magnetised stars should incorporate a boundary condition like ours, but modified to take account of the 
increasing resistivity in the outer regions of the neutron star. 
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APPENDIX A: AXISYMMETRIC MAGNETIC FIELDS 
Al General forms for magnetic field and current 

We wish to see how the assumption of axisymmetry constrains the geometry of the magnetic field and the current; and hence 
also the form of the Lorentz force. Working in cylindrical polar coordinates, we begin with the equilibrium equation for a 
magnetised rotating fluid: 



where we have rewritten ([TJ above by replacing the usual VP/p term with the gradient of the enthalpy H — J dP 1 /p(P') 
and also explicitly written the centrifugal term as the gradient of a scalar. 

If we now take the curl of (|A1|) then by the vector identity V x V/ = (for any scalar field /) we see that 



implying that C/p is also the gradient of some scalar M. Note that VM ■ B = 0, i.e. M is constant along field lines. 




(Al) 




(A2) 
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Next we write B in terms of a streamfunction it, defined through the relations 

D 1 flu _ 1 flu 

Bm — TT" , £>z = — q — 

nj (72: zu Ozu 

— note that these components give a solenoidal magnetic field, V • B = 0, by construction. Hence 
Now comparing the equation with 
we see that B may be written as 



B ro = _l^,^ = l|i (A3) 
zu oz zu ozu 

Lc field, V • B = 0, b 

1 flu „ 1 flu , . 

B = — e ro + + — — e z . (A4) 

zu oz zu ozu 



_ du du . , „. 

Vu x = -— e ro + — - e z , (A5) 
az ozu 



B = — Vu x e + B e . (A6) 
Note that this implies B • Vit = 0, i.e. u is constant along field lines. Recalling that M also has this property, we deduce that 

M = M(u). (A7) 

Next we turn to Ampere's law in axisymmetry: 

47TJ = V x B = -— ^e ro + — — — e + — — (zuB$)e z . (A8) 

flz V ozu I zu ozu 



Now by comparing the poloidal part of the current 
with the quantity 



= "4^£ ( ^ )e - + ^u£u {ujB * )e * (A9) 



we see that 



V(zuB<h) x es = - — (zuB^e^ + ——{vuBAe,, (A10) 
Oz ozu 

jpoi = — — V(xuBj,) x e*. (All) 
iirzu 

Next we consider the toroidal part of the current jtor = j0 e anc ' rewrite using the definition of the streamfunction u: 

dB^ 0B Z 1(3(1 du\ 8 2 u\ 
4tt^ = — — - = w _ ___ +_). A12 

az cro sc7 \ cto \ zu ozu J Oz A J 

For brevity we define a differential operator A* by 

A, = — -j + A13 

ozu z zu ozu oz z 

Now using this definition together with (|A11|) and ()A12[) we see that the current may be written as 

l 1 
47rj = — V(zuB^) x A,u e^. (A14) 

CE7 ZU 

Our two key results from this section so far are the expressions (|A6|) and (|A14I) for the general form of an axisymmetric 
magnetic field and current, respectively. Next we consider the form of the Lorentz force arising from these two quantities. We 
see that in general 

C = j x B = (jpoi + j^e^) x (B poi + B 

= jpoi x Bp 0i + j^e^ x B poi + fi^jpoi x . (A15) 

Returning to our original force balance equation (|A1|I we note that the pressure, gravitational and centrifugal forces are 
axisymmetric (i.e. no ^-dependence); therefore C is also axisymmetric and its toroidal component must vanish: 

Ctor = jpoi x B poi = 0. (A16) 

At this point there are two ways to proceed: either Bp ; is non-zero, in which case B po i and ] po i are parallel; or B po i = 0. We 
shall consider these cases separately in the next two subsections. 



A2 Mixed poloidal and toroidal fields; the Grad-Shafranov equation 

We have shown that the requirement (|A16[) follows from the axisymmetry of our problem. In this subsection we consider the 
case where B po ; and j po ; are parallel, corresponding to a magnetic field with both poloidal and toroidal components. We will 
see that the form of purely poloidal magnetic fields may be found as a special case of the general mixed-field configuration. 
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Recall from (|A"6) and (1A11[) that 
B„ oi = — Vm x e,i 

VJ 

jpoi = -r — V(roB^) x e . 
47rn7 

Knowing that these two quantities are parallel we see that u and vjB^ must be related by some function /: 

vjB^ = /(«). (A17) 

Next we evaluate the non-zero Lorentz force components, i.e. C po i from (|A15[I . Using the pair of equations at the start 
of this subsection, we find that 



and similarly 



x Bp i = x ( — Vu x ] = — (Vu - e^e^ ■ Vm)) = — Vu (A18) 

VJ J VJ VJ 

jpoi x e</> = -— !— VOarBs). (A19) 

47TOT 



Now using these expressions in (|A15[) . together with the relation injj, = — — A*« from (|A14[) . we find that 



■inC = - — A*u X7u - —Bs\7(zuBa) (A20) 

VJ 1 VJ 

which, recalling the definitions VM = Cj p and /(u) = vjB^, becomes 

47rpVM = — —A*u X7u - — f(u)Vf(u). (A21) 

Since M and / are both functions of u alone we are able to rewrite VM(it) and V/(u) using the chain rule, to give 

dM „ 1 . „ 1 „, . df — , , - 

- 47rp— Vu = -^A„ti Vu H t/ u ^Vm. (A22) 



Now provided Vu 7^ we have 



which is the Grad-Shafranov equation i Grad fc Rubin 1958h 



dM 1 /. ,, , d/\ 

= - -J Aji + f(u) -L , (A23) 



We now return to the general form of an axisymmetric current (IA14[I . replacing vjB$ with f(u) and using the chain rule 
to give: 

47rj = — -f-Vu x ej, - — A*ue0. (A24) 
vj du vj 

We now use ()A6[) to make the replacement -Vu x erf, = B po ; and the Grad-Shafranov equation (]A23h to eliminate A*u from 
(fAl4|) : 

47rj = -f B poi + — 47T7i7 2 p— + /(u)-^- • (A25) 



du vj \ du du 

Finally we use the definition / = vjB^ and B = B po ; + B^e^ to yield an expression for the current in terms of the magnetic 
field and the derivatives of the functions M(u) and f(u): 

1 d/„ dM ,. nn . 

J =4^i B + ^^ (A26) 

A3 Purely poloidal field 

Having arrived at an expression for an axisymmetric current associated with a mixed poloidal-toroidal field ()A26[) , we may 
straightforwardly specialise to purely poloidal magnetic fields by choosing f(u) as a constant. Then 4^ = and the mixed 
term vanishes from the expression for j, leaving only a toroidal current 

j = pw—e^ A27 
du 

and hence a purely poloidal field, by Ampere's law. 



A4 Purely toroidal field 

In the previous subsection we showed that (|A26[) may be trivially reduced to the poloidal-field case. However it is clear from 
the form of (|A26[) that there is no choice of / and M which yields a poloidal current (or equivalently a toroidal field). Setting 
M(u) to be a constant, for example, results in the general expression for a force-free field 
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which is of less interest to us, as we aim to study distortions caused by magnetic fields. 

It is clear that the derivation used for mixed fields does not hold in the toroidal-field case. Previously we were able to use 
(|A16I) to simplify the current-field relation, but no such constraint is provided for a toroidal field, where B po ; = 0. Accordingly 
we must return to subsection I A 1 1 where we found that 

Bp i = — Vm x e<k 

zu 

jvol = —^—ViujB^) x 

(from equations (|A6|I and (|A11|) '). Since B po ; = we no longer require vjB^, to be a function of u; indeed the streamfunction 
u will not even enter our final solution. We also recall that the general form of an axisymmetric Lorentz force is given by 
(|A15|) . which in the case of B po ; = reduces to 

£- = B^i po i x e$. (A29) 
Using (|A11[I to replace j po ; in this expression then gives 

1 (VM^)xe^xe^-ivK). (A30) 



4ttvj 4nzu 

Again recalling previous work in this section, we note that taking the curl of (| A 1|) shows that V x (C/p) = 0. We use this 

fact together with the vector identity V x (/V(?) = V/ x Vg to rewrite ()A30[) as 

V (^) x V{vjB^) = 0. (A31) 



If we write j£ in the above expression as —^tziB^ and use the chain rule, some algebra leads to 



50 \7(puj 2 ) x V(roB ) = 0. (A32) 



p 2 CC7 3 



Provided B^/p 2 zu 3 / we then deduce that V(pzu 2 ) x T7(zjBj,) = and hence that pvj 2 and zuB^ are related by some 
function h, i.e. 

vjB^ = h(pzu 2 ). (A33) 

As before we now define a magnetic function M through C/p = VM (note that here M need not be a function of the 
streamfunction u of previous sections). From (|A30|I and ()A33[) we then find that 

h(pza 2 ) , 
Airpuj 2 

By the chain rule we have Vh(y) = where we have introduced the notation 7 = pvj 2 . Given this we have 

VM = -MZ)?V 7 (A35) 
4-7T7 07 



VM = -^^VMP^ 2 ). (A34) 



and so 

rp ™~ hdh 
4-7T J 7 d'y 



M — — - / dry. (A36) 

li" Jo 



